1 Load data

data <- readRDS(list.files(data_path, pattern="FP_matrix", full.names=T))
covs <- readRDS(list.files(data_path, pattern="FP_metadata", full.names=T))
covs <- covs[, c("CERAD", "Braak", "CDR", "plaqueMean")]
protein <- "APP"
index <- match(protein, rownames(data))
input <- data[-index, ]
toPredict <- as.vector(unlist(data[index, ]))
dim(input)

2 Cutoff selection

r <- testAllCutoffs(exprData=input,
                    target=toPredict,
                    covs=covs,
                    train.split=0.7,
                    nfolds=5,
                    t=10,
                    path=paste0(tgcn_path),
                    targetName=protein,
                    tissueName="MSBB",
                    seed=1234,
                    cutoffs=10:1,
                    n=100, 
                    m=10, 
                    s=10, 
                    minCor=0.3,
                    maxTol=3,
                    save=T,
                    overwrite=F)
r <- readRDS(paste0(tgcn_path, "/Net/APP_MSBB_TGCNs_all_cutoffs.rds"))
grid.arrange(r$selectCutoff$nHubs, r$selectCutoff$stats + theme(legend.position="bottom"))

p <- lapply(r$nets, function(cutoff) cutoff$GOenrich$plotStats) 

ggarrange(p$c8 + theme(text=element_text(size=10)) + scale_y_continuous(limits=c(0,25)),
          p$c9 + theme(text=element_text(size=10)) + scale_y_continuous(limits=c(0,25)), 
          p$c10 + theme(text=element_text(size=10)) + scale_y_continuous(limits=c(0,25)), 
          ncol=3, nrow=1, common.legend=T, legend="bottom")

3 TGCN for each cutoff

3.1 Cutoff 8

3.1.1 Modules size selection

r$nets$c8$net$moduleSizeSelectionPlot

3.1.2 Modules correlation

r$nets$c8$net$plotCorr

3.1.3 Modules composition

DT::datatable(r$nets$c8$net$modules)

3.1.4 Cross tab plot

3.1.5 GO enrichment stats

grid.arrange(r$nets$c8$GOenrich$plotStats, r$nets$c8$GOenrich$plotNterms, nrow=2)

DT::datatable(r$nets$c8$GOenrich$terms)

3.1.6 Reduced GO terms at TGCN level

r$nets$c8$GOenrich$plotReduced

reduced_terms <- r$nets$c8$GOenrich$reducedTerms
treemapPlot(reduced_terms)

wordcloudPlot(reduced_terms, min.freq=1, colors="black")

3.1.7 Cell-type enrichment

r$nets$c8$CTenrich$plot

r$nets$c8$CTenrich$sigTests

3.1.8 Module-trait corr

r$nets$c8$moduleTraitCorr$plot_pval

3.2 Cutoff 9

3.2.1 Modules size selection

r$nets$c9$net$moduleSizeSelectionPlot

3.2.2 Modules correlation

r$nets$c9$net$plotCorr

3.2.3 Modules composition

DT::datatable(r$nets$c9$net$modules)

3.2.4 Cross tab plot

3.2.5 GO enrichment stats

grid.arrange(r$nets$c9$GOenrich$plotStats, r$nets$c9$GOenrich$plotNterms, nrow=2)

DT::datatable(r$nets$c9$GOenrich$terms)

3.2.6 Reduced GO terms at TGCN level

r$nets$c9$GOenrich$plotReduced

reduced_terms <- r$nets$c9$GOenrich$reducedTerms
treemapPlot(reduced_terms)

wordcloudPlot(reduced_terms, min.freq=1, colors="black")

3.2.7 Cell-type enrichment

r$nets$c9$CTenrich$plot

r$nets$c9$CTenrich$sigTests

3.2.8 Module-trait corr

r$nets$c9$moduleTraitCorr$plot_pval

3.3 Cutoff 10

3.3.1 Modules size selection

r$nets$c10$net$moduleSizeSelectionPlot

3.3.2 Modules correlation

r$nets$c10$net$plotCorr

3.3.3 Modules composition

DT::datatable(r$nets$c10$net$modules)

3.3.4 Cross tab plot

3.3.5 GO enrichment stats

grid.arrange(r$nets$c10$GOenrich$plotStats, r$nets$c10$GOenrich$plotNterms, nrow=2)

DT::datatable(r$nets$c10$GOenrich$terms)

3.3.6 Reduced GO terms at TGCN level

r$nets$c10$GOenrich$plotReduced

reduced_terms <- r$nets$c10$GOenrich$reducedTerms
treemapPlot(reduced_terms)

wordcloudPlot(reduced_terms, min.freq=1, colors="black")

3.3.7 Cell-type enrichment

r$nets$c9$CTenrich$plot

r$nets$c9$CTenrich$sigTests

3.3.8 Module-trait corr

r$nets$c9$moduleTraitCorr$plot_pval

4 Hub gene module enrichment

go_results <- r$nets$c8$GOenrich$terms
plots <- getReducedTermsPlots(go_results, module=T)
saveRDS(plots, paste0(tgcn_path, "plots_APP_MSBB.rds"))
go_results <- r$nets$c8$GOenrich$terms
sort(table(go_results$query), decreasing=T)
## 
##    UCP2 ATP6AP1  ATP1B1  DHCR24    SDHC    RTN1   LAMP1 HNRNPA0   OTUB2   ABCG1 
##     430      68      64      44      43      35      25       9       8       8 
##   NEK11  NKAIN2   SART1    MCM9    STIL   CAPZB PHYHIPL  ZNF491   RFLNB   RNF32 
##       7       6       5       5       4       4       4       3       2       0 
##    hubs 
##       0
plots <- readRDS(paste0(tgcn_path, "plots_APP_MSBB.rds"))

4.1

plots$ATP1B1$scatterPlot

plots$DHCR24$scatterPlot

plots$SDHC$scatterPlot

plots$UCP2$scatterPlot